args DepVars

clear

global plotopt scheme(s2mono) xlabel(-2(1)2) ylabel(, angle(0)) graphregion(color(gs16) ilc(gs16)) plotregion(lc(gs0)) xtitle("Years since the regulation", size(vlarge)) legend(off) 

foreach depvar in `DepVars' {

	cap gen `depvar' = 0
	run _labels

	local varlab: variable label `depvar'
	estimates use "${direc_estimates}/BJS_`depvar'-dynamic"
	
	* Extract the results
matrix results = e(b)
matrix V = e(V)

* Extract tau and pre coefficients
scalar tau0 = results[1,1]
scalar pre1 = results[1,4]
scalar pre2 = results[1,5]
scalar tau1 = results[1,2]
scalar tau2 = results[1,3]

* Extract standard errors
scalar se_pre1 = sqrt(V[4,4])
scalar se_pre2 = sqrt(V[5,5])
scalar se_tau0 = sqrt(V[1,1])
scalar se_tau1 = sqrt(V[2,2])
scalar se_tau2 = sqrt(V[3,3])

* Calculate 95% confidence intervals
scalar ci_pre1_upper = pre1 + 1.645 * se_pre1
scalar ci_pre1_lower = pre1 - 1.645 * se_pre1
scalar ci_pre2_upper = pre2 + 1.645 * se_pre2
scalar ci_pre2_lower = pre2 - 1.645 * se_pre2
scalar ci_tau0_upper = tau0 + 1.645 * se_tau0
scalar ci_tau0_lower = tau0 - 1.645 * se_tau0
scalar ci_tau1_upper = tau1 + 1.645 * se_tau1
scalar ci_tau1_lower = tau1 - 1.645 * se_tau1
scalar ci_tau2_upper = tau2 + 1.645 * se_tau2
scalar ci_tau2_lower = tau2 - 1.645 * se_tau2

* Use this option to keep pre and post esitimates/CI as they are, and just set tau0=0
scalar pre1_norm = pre1
scalar pre2_norm = pre2
scalar tau1_norm = tau1
scalar tau2_norm = tau2

scalar ci_pre1_upper_norm = ci_pre1_upper
scalar ci_pre1_lower_norm = ci_pre1_lower
scalar ci_pre2_upper_norm = ci_pre2_upper
scalar ci_pre2_lower_norm = ci_pre2_lower
scalar ci_tau0_upper_norm = 0
scalar ci_tau0_lower_norm = 0
scalar ci_tau1_upper_norm = ci_tau1_upper
scalar ci_tau1_lower_norm = ci_tau1_lower
scalar ci_tau2_upper_norm = ci_tau2_upper
scalar ci_tau2_lower_norm = ci_tau2_lower

* Create a temporary file to store the coefficients for plotting
tempfile coeffs
postfile coef str3 time int period float period_effect float ci_lower float ci_upper using `coeffs', replace
post coef ("pre2") (-2) (pre2_norm) (ci_pre2_lower_norm) (ci_pre2_upper_norm)
post coef ("pre1") (-1) (pre1_norm) (ci_pre1_lower_norm) (ci_pre1_upper_norm)
post coef ("tau0") (0) (0) (ci_tau0_lower_norm) (ci_tau0_upper_norm)
post coef ("tau1") (1) (tau1_norm) (ci_tau1_lower_norm) (ci_tau1_upper_norm)
post coef ("tau2") (2) (tau2_norm) (ci_tau2_lower_norm) (ci_tau2_upper_norm)
postclose coef

*local varlab: variable label `depvar'

* Load the coefficients into memory
use `coeffs', clear
replace ci_lower=0 if period==0
replace ci_upper=0 if period==0

twoway rarea ci_lower ci_upper period, astyle(ci) || function y = 0, range(-2 2) lpattern(solid) lcolor(gs8) || scatter period_effect period , xline(0, lpattern(dash)) connect(l) lcolor(black)  mstyle(p1) mcolor(none)   || scatter period_effect period if period <= 0 , msymbol(X) mstyle(p1)  msize(large)  mcolor(black) mlabpos(6) || scatter period_effect period if period > 0 , mstyle(p1) mcolor(black)   mlabpos(12) ||, $plotopt title("`varlab'", size(vlarge)) 

graph export "${direc_output}/BJS_custom_`depvar'.pdf", replace	

}



